      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
     1 RPL,DDSDDT,DRPLDE,DRPLDT,
     2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
     3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
     4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
C
      INCLUDE 'ABA_PARAM.INC'
C
      CHARACTER*80 CMNAME
      DIMENSION STRESS(NTENS),STATEV(NSTATV),
     1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
     3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)
      double precision E, v, lambda, G
      double precision yield0, yieldn, yieldn1,H,W,M,R
      double precision equalPlasticStrain,dEqualPlasticStrain
      INTEGER i,j
      double precision misesEqualStress,misesEqualStressTry
      double precision hydroStaticPressureTry,hydroStaticPressure
      double precision, DIMENSION(NTENS)::deviatoricStressTry,stressTry
      double precision, DIMENSION(NTENS)::deviatoricStress
      INTEGER,save::cn=0
      cn=cn+1

C      user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCD
C      and, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDT
      E=PROPS(1)
      v=PROPS(2)     
      yield0=PROPS(3)

      lambda=E*v/((1+v)*(1-2*v))
      G=E/(2*(1+v))
      H=(PROPS(5)-PROPS(3))/(PROPS(6)-PROPS(4))

C     Calculate the material Jacobian matrix
      do i=1,NTENS
         do j=1,NTENS
            DDSDDE(i,j)=0
         enddo
      enddo
      do i=1,NDI
         do j=1,NDI
            DDSDDE(i,j)=lambda
         enddo
         DDSDDE(i,i)=lambda+2*G
         DDSDDE(i+NDI,i+NDI)=G
      enddo

C     Calculate stress-try  
      do i=1,NTENS
         stressTry(i)=STRESS(i)
      enddo    
      do i=1,NTENS
         do j=1,NTENS
            stressTry(j)=stressTry(j)+DDSDDE(j,i)*DSTRAN(i)
         enddo
      enddo
     
      call GetMisesEqualStress(stressTry,NTENS,NDI,misesEqualStressTry
     1   ,deviatoricStressTry,hydroStaticPressureTry)


C     get  the nth equalPlasticStrain from STATEV
      equalPlasticStrain=STATEV(1)

C     Calculate  the nth yield
      yieldn=yield0+H*(equalPlasticStrain)

         call GetMisesEqualStress(STRESS,NTENS,NDI,misesEqualStress
     1   ,deviatoricStress,hydroStaticPressure)
         if (yieldn.eq.yield0.and.misesEqualStressTry.gt.misesEqualStress ) then
            
            R=(yield0-misesEqualStress)/(misesEqualStressTry-misesEqualStress)
            do i=1,NTENS
               do j=1,NTENS
                  STRESS(i)=STRESS(i)+R*DDSDDE(i,j)*DSTRAN(j)
               enddo
            enddo
            call GetMisesEqualStress(STRESS,NTENS,NDI,misesEqualStress
     1         ,deviatoricStress,hydroStaticPressure)
            do i=1,NTENS
               DSTRAN(i)=(1-R)*DSTRAN(i)
            enddo
         endif
C     Calculate dEqualPlasticStrain

        dEqualPlasticStrain=0
        M=3*G/(misesEqualStress*(H+3*G))
        do i=1,NTENS
            dEqualPlasticStrain=dEqualPlasticStrain+M*deviatoricStress(i)*DSTRAN(i)
        enddo
C       update jacbian matrix        
        w=9*G**2/(misesEqualStress**2*(H+3*G))
        do i=1,NTENS
            do j=1,NTENS
                DDSDDE(i,j)=DDSDDE(i,j)-w*deviatoricStress(i)*deviatoricStress(j)
            enddo
        enddo
C     Calculate stress     
        do i=1,NTENS
            do j=1,NTENS
            STRESS(j)=STRESS(j)+DDSDDE(j,i)*DSTRAN(i)
            enddo
        enddo
         
      else
C     Calculate stress  
        do i=1,NTENS
            do j=1,NTENS
            STRESS(j)=STRESS(j)+DDSDDE(j,i)*DSTRAN(i)
            enddo
        enddo        
         
        dEqualPlasticStrain=0
      end if

      do i=1,NTENS
        STRAN(i)=STRAN(i)+DSTRAN(i)
      enddo
    

      STATEV(1)=equalPlasticStrain+dEqualPlasticStrain
      RETURN
      END
      SUBROUTINE GetMisesEqualStress(STRESS,NTENS,NDI,misesEqualStress
     1   ,deviatoricStress,hydroStaticPressure)
      INCLUDE 'ABA_PARAM.INC'
      DIMENSION STRESS(NTENS)
      double precision, DIMENSION(NTENS)::deviatoricStress
      double precision misesEqualStress,hydroStaticPressure
      INTEGER NTENS,NDI
      hydroStaticPressure=0
      do i=1,NDI
         hydroStaticPressure=hydroStaticPressure+STRESS(i)
      enddo
      hydroStaticPressure=hydroStaticPressure/3
      do i=1,NDI
         deviatoricStress(i)=STRESS(i)-hydroStaticPressure
      enddo
      do i=1+NDI,NTENS
         deviatoricStress(i)=STRESS(i)
      enddo
      misesEqualStress=0
      do i=1,NDI
         misesEqualStress=misesEqualStress+deviatoricStress(i)**2
      enddo
      do i=1+NDI,NTENS
         misesEqualStress=misesEqualStress+2*deviatoricStress(i)**2
      enddo
      misesEqualStress=sqrt(1.5*misesEqualStress)
      END      